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We extend Kubo's Linear Response Theory (LRT) to periodic input signals with arbitrary shapes 
and obtain exact analytical formulas for the energy dissipated by the system for a variety of signals. 
These include the square and sawtooth waves, or pulsed signals such as the rectangular, sine and im- 
pulses. It is shown that for a given input energy, the dissipation may be substantially augmented by 
exploiting different signal shapes. We also apply our results in the context of magnetic hyperthermia, 
where small magnetic particles are used as local heating centers in oncological treatments. 



Kubo's linear response theory [TJ (LRT) has seen an 
immense success in the past decades. It has served as 
the starting point for much of the advances in statistical 
physics and stochastic process, while also acting as the 
foundation to a variety of applications [H [3] . One of the 
many important breakthroughs it provides is the ability 
to relate the dynamic properties of the system with it's 
quasi-static response. This is of particular importance in 
dielectric and magnetic media [21 [3] , where the response 
to a harmonic electric or magnetic field may be accurately 
related to the relaxation time of the system, provided the 
field amplitude is sufficiently small. Relevant examples 
include the electric response of polar molecules and liquid 
crystals [3] , or the magnetic response of small magnetic 
particles [21 0|. The latter, in particular, has seen re- 
newed interest in recent years due to it's potential appli- 
cation in oncological treatments via a technique known 
as magnetic hyperthermia [51 [BJ. In it, with the purpose 
of thermally lysing tumorous cells, one exploits the heat 
dissipated by magnetic nanoparticles under the influence 
of an external high frequency magnetic field. For a har- 
monic stimuli, the dissipation is directly related to the 
imaginary part of the complex susceptibility. Thus, one 
may employ the LRT as a tool in the quest to optimize 
the heat dissipated by the particles. 

The LRT, however, is not restricted to harmonic in- 
puts and may readily be extend to describe the system's 
response to any periodic stimuli expandable in a Fourier 
series [7]. Albeit straightforward, to our knowledge such 
development has not yet been performed analytically. 
Clearly, it is possible to envisage several situations where 
these results may be of use. We will, however, focus pri- 
marily on magnetic hyperthermia. The main motivation 
is that, in maximizing the dissipation of the particles, it is 
also necessary to maintain a sufficiently low exciting fre- 
quency in order to avoid the formation of eddy currents 
inside the patient's body [BJ. Thus, different alternatives 
must be exploited, the field's shape being one possibility. 

It is the purpose of this paper to extend the LRT to 
encompass arbitrary periodic signals. First, a general for- 
mula to numerically evaluate the response is developed. 
Thereafter, we focus primarily on the energy A dissi- 
pated by the system. For a variety of signals [cf . Fig. [Tj , 
we compute exact analytical formulas for A in terms of 



x = lot, where w is the exciting frequency and t is the re- 
laxation time of the system. We begin by studying simple 
signals such as the square and sawtooth waves. Subse- 
quently, we turn to periodic pulsed signals and develop 
formulas for the rectangular and sine pulses. From these, 
the J-pulse response is also readily obtained. Finally, we 
apply these results in the context of small magnetic par- 
ticles and compare them to numerical simulations of the 
magnetic Langevin equation [2j. 

The computations that follow are valid for any in- 
put/output system to which the LRT applies. For con- 
creteness, however, we shall focus on the response of a 
magnetic system with uniform magnetization M(t), to a 
periodic magnetic field H{t) decomposable in a Fourier 
series. We write H(t) = H rj{t) where, 

oo oo 

r/(t) = ciq + (a n cos ncut + b n sin nojtj = c n e lult 

n—l n— — oo 

(1) 

For comparative purposes we normalize the energy in rj (t) 
to match that of a harmonic input [r](t) — coswt]; i.e., 
we write 

ui J \v(t)\ 2 dt = 7T (2) 

period 

Thus, the input energy is now always equal to ttHq and 
the criteria for the applicability of the LRT becomes re- 
stricted exclusively to Hq. 

Assuming that this criteria is satisfied, we expect that 
the magnetization will respond linearly to H(t), albeit 
possibly with a phase lag. Whence, we write 

t 

M(t)= f H(t') X (t-t')dt, (3) 



where x(t) is a function describing the system's response. 
Inserting Eq. in Eq. @ we obtain 

oo 

M(t) = H c n x(noj)e m ^, (4) 
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where 



X(nw) 



e- in " t X (t)dt = xa 



C{t)e 



dt 



(5). 



The second equality follows directly from the Kubo re- 
lation pp. Here xo — x(0) is the static susceptibility 
and C(t) is the autocorrelation function of the system. 
Knowledge of these two quantities enables one to employ 
Eqs. 0, Q and ([5| to compute the response of the 
system to any periodic stimuli. 

A parametric plot of (H(t),M(t)) yields a hysteresis 
loop, the area of which is precisely the average energy 
dissipated per cycle (which follows from the first law of 
thermodynamics): A = jHdAI. Let us write dM = 
M(t) dt. Then, since M(t) is a continuous function of 
time [H(t) need not be], we may differentiate it's Fourier 
series term by term [7j. Thus, M(t) will also be described 
by a Fourier series and we may use Parseval's theorem 
[7] to write 



(6) 



71=1 



This formula is entirely general, valid for any autocorre- 
lation function. 

Henceforth, we specialize the calculations further to 
systems described by a Fokker-Planck equation. In this 
case C (t) is described by an infinite sum of decaying ex- 
ponentials, each representing a possible relaxation mech- 
anism of the system. Often, however, a single exponential 
provides the dominant contribution; i.e., C(t) = e~*/ T , 
where r is then referred to as the relaxation time of the 
system. In what follows we consider only such form for 
C(t). This assumption, however, is not restrictive given 
the linearity of all the equations involved: if C(t) is de- 
scribed by more than one exponential, their contributions 
may simply be appended to Eq. |6]) , weighted with proper 
coefficients. For C(t) = e - */ r , Eq. (|5| becomes 



X(nw) = xo- 



(7) 



Note that, now, xi nuJ ) is the only quantity where the 
frequency has a definite influence. Everywhere else we 
may equivalently set oj = 1, making the fundamental 
period of the signal equal to 2ir; thus, we henceforth take 

t G [-7T,7r]. 

Since the average input energy per cycle is ttHq [cf. 
Eq. we may define the efficiency in converting elec- 
tromagnetic energy into thermal energy — or, more gen- 
erally, the input /output energy gain — as = A/itHq. 
Whence, inserting Eq. ^ into Eq. ^ we obtain: 
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n x 



n=l 



1 + (nx) 2 



LOT 



(8) 



This formula gives the dissipation efficiency for any sys- 
tem described by a single relaxation time. As mentioned, 
the generalization for more than one relaxation time is 
straightforward. It is quite remarkable that, except for 
XOj all other properties of the system condense into a sin- 
gle variable: x = lot. For simplicity, we henceforth set 
Xo = 1- 

As expected, for a harmonic signal [r)(t] — cost], 
Eq. ([8]| yields: 



Q(x) = 



1 



(harmonic wave) 



(9) 



This result is plotted in Fig. [2j curve 1, where it is seen 
to have a maxima at x = lot = 1 . 

One point requires further clarification: suppose we 
choose r](t) in Eq. ^ such that = 1 and all other co- 
efficients are zero (i.e., rj(t) = cos(2t)). This would then 
yield a result which is twice that obtained by replacing 
x with 2x in Eq. (I9j), a consequence of the n 2 term in 
Eq. The reason for this apparent contradiction is 

that Q(x) describes the energy conversion efficiency per 
period for a signal whose fundamental period is 2ir (or 
2tt/lo in real units). Hence the factor of two in this ex- 
ample would follow from counting the energy dissipated 
in two fundamental periods instead of one. A very im- 
portant consequence follows from this argument: take, 
for concreteness, the point x = 1 corresponding to the 
maxima of Eq. |9). Since the corresponding maxima of 
higher order harmonics occur at different values of x, 
one may argue that replacing the main harmonic with 
a weighted sum containing higher harmonics should al- 
ways reduce the net dissipation. However, albeit having 
a smaller dissipation per fundamental period, the higher 
harmonics are oscillating with respect to the fundamen- 
tal period of the first harmonic. In other words, the n-th 
order harmonic completes a total of n periods during a 
time interval of 2ir, thence compensating for it's lower 
dissipation. 

In Fig. [T] we summarize the signals investigated in this 
paper. We begin with the square wave [Fig. [l] (a)]: 



v(t) = 



1 



1*1 <f 
1 otherwise 



(10) 



where the factor 1 / \/2 was required in order to satisfy 
Eq. ([2]). Inserting the corresponding Fourier coefficients 
in Eq. (j8J), and using the partial fraction expansion of 
tanh(z)/z, we obtain 



2 / 7T 

Vt(x) = — tanh ( — 
7T \2x 



(square wave) 



(11) 



This result is plotted in Fig. [5J curve 2. As can be seen, 
it is more efficient than the harmonic signal for all values 
of x (i.e., it dissipates more energy for a given input en- 
ergy). It also presents a "plateau" below x = 1. This is 
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important given that real systems always have some dis- 
tribution of relaxation times. Thus, with a square wave 
the dissipation is expected to be more homogeneous. The 
fact that fl(x — > 0) — > 2/ir, which is clearly unphysical, 
is a consequence of the discontinuity of the input signal; 
evidently, in real systems Cl(x — > 0) — > since the input 
is produced by a finite number of harmonics. To further 
emphasize this point we present on the inset of Figj2]thc 
result of numerically evaluating the sum [cf. Eq. ([8|] up 
to 10 (open circles) and 100 (filled circles) harmonics. As 
can be seen, increasing the number of harmonics signifi- 
cantly enhances the tendency of the curve to remain flat 
close to x = 0. 

Next, we turn to the sawtooth wave [Fig. [I] (b)]: 



(12) 



Inserting the corresponding Fourier coefficients in Eq. ([8]) 
and using, this time, the partial fractions expansion of 
coth(z)/z, we obtain 

3 r 



coth 



(sawtooth wave) (13) 



This result is plotted in Fig. [5] curve 3. It dissipates 
more than both the harmonic and square waves, and has 
the zero-frequency limit ft(x — >• 0) — > 3/tt, again due to 
the discontinuity. 

An important aspect of Eq. Q is that the signal's 
shape enters only in terms of the combination (a^ + 6^). 
This means that entirely different signals may have the 
exact same efficiency. An example is the signal in 
Fig. [I] (c), corresponding to r](t) = log [2 cos(i/2)]. It's 
dissipative properties coincide exactly with those of the 
sawtooth wave. 

A closed form solution for the (continuous) triangular 
wave also exists. Such, however, is of little practical inter- 
est since it nearly coincides with the harmonic efficiency 
[cf. Eq. j9} and Fig. [2| curve 1]. This is expected given 
the extremely rapid convergence of it's Fourier series. 

The step-ladder signal depicted in Fig. [I] (d) approx- 
imately mimics real field variations. It can be written 
as 



1 



a 



/c=i 



u \t + 



kit 



a 



u t 



a ) 



,(14) 



where u is the unit-step function and a = 2, 3, 4, . . . de- 
fines the number of steps, with 2 referring to the square- 
wave. The signal in Fig. [I] (d) is for a — 4. In the limit 
a — > oo we recove r the harmonic field. A normalization 
constant y^ 2(a+i) is also missing. No closed form solu- 
tion exists for this signal. Notwithstanding, the corre- 
sponding sum may always be evaluated numerically. Re- 
sults for a = 4 are shown in Fig.[2j curve 4. As expected, 
it lies between the square and harmonic waves, illustrat- 
ing a gradual transition taking place between these two 
asymptotes. 



Next we turn to pulsed signals. We begin with the 
rectangular pulse [Fig{T](e)]: 



rj(t) = 




1 Itl < z 

ii — a 

otherwise 



(15) 



Here a represent the width of the pulse: a = 1 corre- 
spond to a straight line and a — > oo to a 5-pulse. The 
function in Fig. [I] (e) is for a = 4 and the square-wave 
[Eq. ( 10 )] correspond to a = 2. However, the signal is 



no longer symmetric with respect to r\ = 0, which means 
that energy is being wasted in the do term of the Fourier 
series in Eq. ([!]). It thus follows that the efficiency for 



a = 2 will be half of that given by Eq. ( 11 ). 



The simplest way to obtain a general formula for the 
rectangular pulse [Eq. (15)], valid for all a > 1, is by 
induction. For instance, when a = 4 the result is H.4 (x) = 
(l/7r)[tanh(7r/2:c) + tanh(7r/4a;)]. The similarity between 
this result and Eq. ( |TT| ) incites the idea that the induction 
formula may be written as a sum of hyperbolic tangents. 
Unfortunately, this is not the case. However, if rewrite 
S^2 and f24 in the form 



2 sinh(7r/2x) sinh(7r/2x) 
7r sinh(7r/a;) 

4 sinh(7r/4a;) sinh(37r/4:r] 
7r sinh(7r/x) 



then another pattern becomes clearly visible inciting us 
to write: 



asinh(iz) sinWs^in 

Q a (x) = Kax! , : a XJ 

y ' 7T sinh(z) 



(16) 



(rectangular pulse). It turns out that this result is ac- 
tually valid for all values of a (with a > 1), including 
non-integers. Clearly, the correctness of the formula is 
easily tested by comparing it with the numerical calcula- 
tion of the sum in Eq. Q. 

Results for a = 2 (square wave), 4, 6 and 8 are shown 
in Fig. [3] together with the harmonic response, shown in 
dashed for comparison. As can be seen, the narrower 
the pulse (larger a), the more efficient is the dissipation. 
For large x the function behaves as (1 — (l/a))/x, which 
is smaller than the harmonic efficiency (which goes as 
1/x). At the other extreme, close to x = 0, we have 
that £l a (x — » 0) — » a/2n; i.e., it scales linearly with a. 
Taking the limit a — > 00 we obtain the (5-pulse response 
^oo(^) = 1/x, which is illustrated in Fig. [3]in a dotted 
line. To obtain the efficiency for a pulse symmetric with 
respect to 77 = 0, one need only multiply Eq. ( fl6] ) by 2. 
In this case it is worth noting that for all pulses with 
a > 2, fl a remains above the harmonic efficiency for all 
x. 

Next we turn to the more realistic sine pulse [Fig.[l](f )] : 



v(t) = 



cos(crf/2) |i|<s 







otherwise 



(17) 
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The function in Fig. [ljf) is for a = 4. We restrict 
our analysis to a > 1, corresponding to the pulses with 
t?(±7t) = 0. The calculations for this signal are somewhat 
more cumbersome. As before, the simplest approach is 
by induction. For conciseness, we give only the result: 



Q a (x) = 



[4 + (ax) 2 



4 + (ax) 2 



(18) 



8ax cosh (g^f) cosh 



a x / 



sinh ( ^ 



(sine pulse) . Similarly to Eq. ([16]) , this formula turns out 
to be valid for all a > 1. 

Results for a = 2, 4, 6 and 8 are shown in Fig.|4j whose 
details are similar to Fig. [3] As can be seen, it shares a 
clear similarity with the rectangular pulse depicted in 
the latter. However, since it is continuous, we now have 
Sl(x — > 0) — > 0. As before, sufficiently narrow pulses 
(a > 4) are considerably more efficient than the harmonic 
wave when x < 1. 

We have thus far focused on an autocorrelation func- 
tion of the form C(t) = e~ l l T . Near a critical point, how- 
ever, the decay is known to become algebraic: C(t) oc 
t~ s . For 8 < 1 we may use Eq. (pM) to show that 
ft oc u! 5 for all signal shapes; i.e., the only thing chang- 
ing from one signal to another is the proportionally con- 
stant. Slight care must be taken, however, in the fact 
that the Fourier sum now reads (a£ + l%)n s+1 and may 
thus present convergence issues. 

Finally, we apply these results to the problem of mag- 
netic hyperthermia 5 , 6]. The magnetic Langevin equa- 
tion and it's corresponding Fokker-Planck equation were 
first introduced in Ref. 4] and are described in detail in 
Ref.[2J Conveniently, over a broad frequency interval (the 
ferromagnetic resonance region excerpted), the autocor- 
relation function is adequately described by a single de- 
caying exponential with relaxation time r(a) — ^ \p^e° ■ 
Here, To ~ 1CP 9 s and a is the ratio between the en- 
ergy barrier separating the stable energy minima and 
the thermal energy H]: a = Kv/ksT, where K is the 
anisotropy constant and v is the particle's volume. On 
the other hand, the static susceptibility may be written 
as xo — const x (a — 1). Finally, we have also fixed 
ujtq = 1CP 4 for concreteness; the response to other fre- 
quencies is qualitatively similar. 

In Fig. [5] we present in green lines (filled circles) the 
efficiency as a function of a for the harmonic, square and 
sawtooth waves (from inner to outermost). We also show 
in red (empty circles) the calculations for the sine pulse 
[Eq. ( 18 )] with a — 2, 4 and 8 (again, from inner to outer- 
most). Graphing ft vs. a enable us to directly related the 
efficiency to the fundamental magnetic parameter of the 
system. It can be seen that there is a maxima associated 
with each signal, corresponding to the value of a to which 
the particles should be tailored in order to maximize the 



dissipation. As for the asymptotic behavior of f2, when 
a is large all functions tend to zero exponentially. For 
small a, on the other hand, the continuous signals (har- 
monic and sine pulse) scale roughly as f2 oc <r 4 whereas 
the discontinuous square and sawtooth waves present a 
linear dependence, oc a. 

The scattered points in Fig. [5] were computed di- 
rectly from the magnetic Langevin equation, first ob- 
taining M(t) from the numerical solution of a hierar- 
chy of differential recurrence relations, and then the area 
[A = fHdM] by numerical integration. The computa- 
tional details are described thoroughly in Ref. [H As can 
be seen, the agreement between both methods — which 
are of entirely different nature — clearly establish the 
correctness of the formulas presented in this paper. 

Finally, we acknowledge the experimental challenges 
of producing non-harmonic AC magnetic fields for hy- 
perthermia. The technique presently used is based on 
LC resonant circuits; whence, one alternative would be 
to stack synchronized circuits to produce the necessary 
harmonics. However, other alternatives may also exist for 
particular signals, such as the sawtooth wave which re- 
quires a ramp-like field increase, or general pulsed fields 
created from current pulses. Another important topic 
regards the possible biological side effects of employing 
higher order harmonics (related to the formation of eddy 
currents). Unfortunately, we are at present unable to 
provide an adequate answer to this question in view of 
the lack of experimental results on the subject. We note, 
however, that the formation of eddy currents is propor- 
tional to the product ujHq. For this reason, the effect of 
higher order harmonics is expected to be, at least par- 
tially, mitigated by their smaller amplitudes (the Fourier 
coefficients always decay faster than 1/n). 

In conclusion, we have shown how Kubo's linear re- 
sponse theory may be extended to account for arbitrary 
signal shapes. The energy dissipated by the system was 
studied for several common signals and analytical formu- 
las were obtained in a variety of cases. It was shown that 
for the same energy input, substantial improvements in 
the dissipated output can be realized using different sig- 
nal shapes. Even though the development was performed 
with respect to a magnetic system, the calculations here 
presented are entirely general and are expected to remain 
valid for any system where the linear response theory is 
applicable. By comparing the exact calculations with nu- 
merical simulations of magnetic hyperthermia we have (i) 
confirmed the exactness of our results and (ii) illustrated 
an important application of non-harmonic stimuli. 
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FIG. 1. Collection of signals studied, (a) square wave, 
Eq. {10}; (b) sawtooth wave, Eq. Q; (c) n(t) = 
log[ 2cos (f/2)]; (d) step-ladder, Eq. ( |14[ |; (e) rectangular pulse, 
Eq. (pf ; (f) sine pulse, Eq. 0\. 
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FIG. 2. Efficiency vs. x for different signals, (curve 1) Har- 
monic wave, Eq. (c urv e 2) Square wave, Eq. (Ill; (curve 
3) Sawtooth wave, Eq. (131; (curve 4) Step-ladder with a = 4 

(pf and fnj]. (inset) Nu- 
for the square- wave 



[evaluated numerically from Eqs. 
merical evaluation of the sum in Eq. 
extending up to 10 (open circles) and 100 (filled circles) har- 
monics. For comparative purposes, solid lines denote the har- 
monic and square- waves, as in the main plot. 
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FIG. 3. Efficiency vs. x for the rectangular pulse [Eq. ( 16 1] for 
different values of a. Harmonic efficiency is shown in dashed 
and 5- impulse response, SI = 1/x, is shown in dotted. 
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FIG. 4. Efficiency vs. as for the since pulse [Eq. ( 18 1] for 
different values of a. Harmonic efficiency is shown in dashed 
and (5-impulse response, Q = 1/x, is shown in dotted. 




10 12 14 



FIG. 5. Efficiency vs. a, the ratio of the energy barrier sepa- 
rating stable energy minima to the thermal energy (see text 
for details), (green lines, open circles) Harmonic, square and 
sawtooth waves (from inner to outermost), (red lines, filled 
circles) Sine pulse with a = 2, 4 and 8 (from inner to out- 
ermost). Scattered points were computed from the Magnetic 
Langevin equation, simulating the magnetization M(t) and 
numerically evaluating the area of the hysteresis loop. 



